########################################################## #
###   R code to analyze data for
###   'When should the majority rule? Experimental
###     evidence for Madisonian judgments in five cultures'
###   by Bor, Mazepus, Bokemper and DeScioli
########################################################## #
###   This code will produce all models, tables & figures in paper & appendix

###   Code by Alexander Bor (alexander.bor@ps.au.dk)
########################################################## #
###   Code written with R version 3.6.2


### SETUP ###

# Install/update packages (if necessary)
# install.packages("lme4")
# install.packages("tidyverse")
# install.packages("pander")
# install.packages("stargazer")
# install.packages("pwr")
# install.packages("here")
# install.packages("brms")
# install.packages("arm")

# Load packages
library(lme4)
library(pander)
library(stargazer)
library(here)
library(pwr)
library(brms)
library(arm)
library(tidyverse)

# custom functions
se <- function(x) {
  xnn <- x[!is.na(x)]  # x with no na
  (sd(xnn) / sqrt(length(xnn)))
}


# make sure import script is run first
# source(here("B. command files/import.R"))

# load pooled, cleaned, long-format data
load(here("C. analysis data/minority_data.Rdata"))

# load data without exclusions
load(here("C. analysis data/minority_all_data.Rdata"))

######################### #
# Sample characteristics  #################################################################
######################### #

# This code reproduces sample sizes reported in main text and
#   and Table A1 in the online appendix
sample_char <- df %>%
  # we have multiple obs per resp. drop these.
  filter(!duplicated(ID)) %>%
  group_by(country) %>%
  # summarise sample N, gender ratio, age and ideology
  summarise(
    n = n(),
    sex = sum(sex == 2, na.rm = T) / n,
    age.m = mean(age, na.rm = T),
    age.sd = sd(age, na.rm = T),
    ideo.m = mean(ideo, na.rm = T),
    ideo.sd = sd(ideo, na.rm = T)
  )

# add info on sampling method
sample_char$sampling <-
  c(
    "University students",
    "University students",
    "Mechanical Turk",
    "University students",
    "Mechanical Turk"
  )

# reshuffle vars
sample_char <- sample_char %>%
  dplyr::select(country, sampling, n, sex, age.m, age.sd,
                ideo.m, ideo.sd)

panderOptions("digits", 3)
pander(sample_char)

write.csv(sample_char, file = here("D. documents", "appendix_table_a1.csv"))


# Statistics reported in Method / Participants
sum(sample_char$n) # number of participants in sample
nrow(df) # number of observations

# Participants excluded for failing comprehension check
length(unique(df_all$ID)) - length(unique(df$ID))

(length(unique(df_all$ID)) - length(unique(df$ID))) / length(unique(df_all$ID))


######################### #
# Descriptive stats of choice ################################################################
######################### #


# ~ calculate the proportion of choices by country and condition ##############

df_prop_c <- df %>%
  filter(!is.na(decision)) %>%  # exclude NAs from decision variable
  mutate(condition = fct_recode(condition, Control = "control",
                                Minority = "treatment")) %>%
  group_by(condition, country) %>%
  summarize(
    vote = sum(decision == "vote"),
    consensus = sum(decision == "consensus"),
    leader = sum(decision == "leader"),
    chance = sum(decision == "chance"),
    trials = n()
  ) %>%
  gather(decision, frequency, vote:chance) %>%
  mutate(
    decision = fct_relevel(decision, "consensus", "leader", "chance"),
    country = fct_recode(
      country,
      DNK = "Denmark",
      HUN = "Hungary",
      IND = "India",
      RUS = "Russia"
    ),
    p = frequency / trials
  )

# Baseline preference for voting
# Values from this table are quoted in the second paragraph of the results section
baseline_vote <-
  filter(df_prop_c, condition == "Control", decision == "vote") %>%
  arrange(p) %>%
  dplyr::select(condition, decision, country, frequency, trials, p)
pander(baseline_vote)

# determine each participant’s percentage of choosing each decision rule
df_pool <- df %>%
  # group by respondents
  group_by(ID) %>%
  # average outcomes of scenarios
  mutate(
    vote = mean(decision == "vote", na.rm = T),
    consensus = mean(decision == "consensus", na.rm = T),
    leader = mean(decision == "leader", na.rm = T),
    chance = mean(decision == "chance", na.rm = T),
    app_vote = mean(app_vote, na.rm = T),
    app_consensus = mean(app_consensus, na.rm = T),
    app_leader = mean(app_leader, na.rm = T),
    app_chance = mean(app_chance, na.rm = T)
  ) %>%
  # drop duplicated responses for each individual
  filter(!duplicated(ID))

# ~ anova for baseline level of vote by country
summary(aov(vote ~ country, df_pool[df_pool$condition == "control", ]))


# ~ Plot Figure 1. ##################
# This code reproduces Figure 1 in the main text.
ggplot(df_prop_c, aes(
  x = fct_rev(condition),
  y = p,
  fill = decision
)) +
  geom_col() +
  facet_grid(country ~ .) +
  # manual color codes help with unconventional order
  scale_fill_manual(values = c("#606060", "#bcbaba", "#969696", "#353535")) +
  scale_y_continuous(breaks = seq(0, 1, 0.25), labels = scales::percent) +
  coord_flip() +
  xlab("") +
  ylab("") +
  theme_minimal() +
  guides(fill = guide_legend(reverse = T)) +
  theme(
    legend.title = element_blank(),
    legend.position = "bottom",
    legend.spacing.x = unit(0.6, 'cm'),
    panel.grid.major.y = element_blank(),
    text = element_text(size = 13)
  )
ggsave(here("D. documents/fig_1.jpg"), width = 6.5, height = 4.2)


########################## #
# Baseline level of vote by country ###########################################
########################## #

# Values from this table are quoted in the 2nd paragraph of the results section
# They show the distribution of preferences in control vs treatment conditions

df_pool <- df %>%
  group_by(ID) %>%
  mutate(
    vote = mean(decision == "vote", na.rm = T),
    consensus = mean(decision == "consensus", na.rm = T),
    leader = mean(decision == "leader", na.rm = T),
    chance = mean(decision == "chance", na.rm = T),
    app_vote = mean(app_vote, na.rm = T),
    app_consensus = mean(app_consensus, na.rm = T),
    app_leader = mean(app_leader, na.rm = T),
    app_chance = mean(app_chance, na.rm = T)
  ) %>%
  filter(!duplicated(ID))

# ANOVA for variability across countries
summary(aov(vote ~ country, df_pool[df_pool$condition == "control", ]))

########################## #
# Average treatment effect on vote in pooled sample ##########################
########################## #


t.test(vote ~ condition, df_pool)

# ~ country level treatment effects #######

# formula to export stats of interest from t-tests
country.t <- function(data, dv) {
  t <- t.test(formula(substitute(dv ~ condition)), data)
  data.frame(
    meandiff = round(diff(t$estimate), 2),
    t = t$statistic,
    df = t$parameter,
    p = t$p.value
  )
}

by(df_pool, df_pool$country, function(x)
  country.t(x, vote))

# ~ treatments effect conditional on country ######
summary(aov(vote ~ country * condition, df_pool))


########################## #
# Average treatment effect on consensus in pooled sample ######################
########################## #

t.test(consensus ~ condition, df_pool)

# ~ country level treatment effects #######
by(df_pool, df_pool$country, function(x)
  country.t(x, consensus))

# ~ treatments effect conditional on country ######
summary(aov(consensus ~ country * condition, df_pool))

########################## #
# Appendix A3 - Power calculations ##############################################
########################## #

# function to calculate chisq.test effect size "w"
#       from test statistic and number of observations
# https://ncss-wpengine.netdna-ssl.com/wp-content/themes/ncss/pdf/Procedures/NCSS/Chi-Square_Effect_Size_Calculator.pdf
w <- function(chisq, n) {
  sqrt(chisq / n)
}


# Calculate required sample size to get 80% power for forced choice variables
# inputs are estimates reported in DeScioli and Bokemper’s (2018) Study 3
pwr.chisq.test(w(39.4, 117),
               df = 3,
               sig.level = 0.05,
               power = 0.8)
pwr.chisq.test(w(12.13, 117),
               df = 3,
               sig.level = 0.05,
               power = 0.8)
pwr.chisq.test(w(6.82, 117),
               df = 3,
               sig.level = 0.05,
               power = 0.8)

# function to calculate Cohen's d from reported stats (t stat & df)
d <- function(t, df) {
  2 * t / sqrt(df)
}

# Calculate required sample size to get 80% power for appropriateness variables
# inputs are estimates reported in DeScioli and Bokemper’s (2018) Study 3
pwr.t.test(d = d(7.71, 115),
           power = 0.8,
           sig.level = 0.05)
pwr.t.test(d = d(3.98, 115),
           power = 0.8,
           sig.level = 0.05)
pwr.t.test(d = d(4.03, 115),
           power = 0.8,
           sig.level = 0.05)

########################## #
# Appendix C – Additional analyses of participants’ choice of decision rule #########
########################## #

# Plot preference for voting in control group by country against
#       Freedom House democracy index

# create dataframe with FH democracy scores
demscore <- data.frame(
  country = c("DNK", "HUN", "IND", "RUS", "USA"),
  country_full = c("Denmark", "Hungary", "India", "Russia", "USA"),
  demscore = c(97, 76, 77, 20, 89)
)

# ~ Figure C1.############
df_prop_c %>%
  filter(decision == "vote" & condition == "Control") %>%
  right_join(demscore) %>%
  ggplot(aes(x = demscore, y = p)) +
  geom_smooth(method = "lm", se = F, color = "black") +
  geom_label(aes(label = country)) +
  scale_y_continuous(labels = scales::percent) +
  xlab("Freedom House Democracy Score (2017)") +
  ylab("Proportion choosing voting")                    +
  theme_bw()
ggsave(here("D. documents", "fig_c1.jpg"), width = 6.5)

# calculate correlation.
corrdata <- df_prop_c %>%
  filter(decision == "vote" & condition == "Control") %>%
  right_join(demscore)

cor.test(corrdata$demscore, corrdata$p)

# ~ Figure C2. #################################

# calculate the proportion of each answer across conditions and countries
df_prop <- df %>%
  filter(!is.na(decision)) %>%  # exclude NAs from decision variable
  mutate(
    condition = fct_recode(condition, Control = "control",
                           Minority = "treatment"),
    scenario = fct_recode(
      scenario,
      Dinner = "dinner",
      Activity = "activity",
      Company = "company"
    )
  ) %>%
  group_by(condition, scenario, country) %>%
  summarize(
    vote = sum(decision == "vote"),
    consensus = sum(decision == "consensus"),
    leader = sum(decision == "leader"),
    chance = sum(decision == "chance"),
    trials = n()
  ) %>%
  gather(decision, frequency, 4:7) %>%
  mutate(
    decision = fct_relevel(decision, "consensus", "leader", "chance"),
    country = fct_recode(
      country,
      DNK = "Denmark",
      HUN = "Hungary",
      IND = "India",
      RUS = "Russia"
    ),
    p = frequency / trials
  )


ggplot(df_prop, aes(
  x = fct_rev(condition),
  y = p,
  fill = decision
)) +
  geom_col() +
  facet_grid(country ~ scenario) +
  # manual color codes help with unconventional order
  scale_fill_manual(values = c("#606060", "#bcbaba", "#969696", "#353535")) +
  scale_y_continuous(breaks = seq(0, 1, 0.5), labels = scales::percent) +
  coord_flip() +
  xlab("") +
  ylab("") +
  theme_minimal() +
  guides(fill = guide_legend(reverse = T)) +
  theme(
    legend.title = element_blank(),
    legend.position = "bottom",
    legend.text = element_text(size = 10),
    panel.grid.major.y = element_blank(),
    axis.text.x = element_text(size = 9),
    panel.spacing = unit(0.75, "lines")
  )

ggsave(here("D. documents/fig_c2.jpg"), height = 3.8, width = 6.5)


# ~ Multinomial logistic regression model #########################
# run baseline model only with random intercepts
# warning: this may take several minutes
multi0 <- brm(
  decision ~ condition +
    (1 | ID) +
    (1 | country_scenario),
  data = df,
  family = "categorical",
  cores = 3,
  prior = c(set_prior("normal (0, 8)"))
)

# Save model outputs to save time
save(multi0, file =  here("C. analysis data/multilevel_multinom0.Rdata"))
# load("C. analysis data/multilevel_multinom0.Rdata")

# Run main model with random intercepts + random slope for condition
# warning: this may take several minutes
multi1 <- brm(
  decision ~ condition +
    (1 | ID) +
    (1 + condition | country_scenario),
  data = df,
  family = "categorical",
  cores = 3,
  prior = c(set_prior ("normal (0, 8)"))
)

# Save model outputs to save time
save(multi1, file = here("C. analysis data/multilevel_multinom1.Rdata"))
load(here("C. analysis data/multilevel_multinom1.Rdata"))

# plot(multi1) # diagnostic plots. Make sure MCMC chains converge on same value. looks good!
summary(multi1)

# test if added complexity improves model fit with leave-one-out cross-validation.
loo(multi0, multi1) # lower value = better fit!


# Export population level effects.
# Table C1 in Online appendix
pop_results <- data.frame(summary(multi1)$fixed)
names(pop_results)[3:4] <- c("Lower", "Upper")
rownames(pop_results) <-
  c(
    "Consensus - Intercept",
    "Leader - Intercept",
    "Chance - Intercept",
    "Consensus - Vulnerable minority",
    "Leader - Vulnerable minority",
    "Chance - Vulnerable minority"
  )

pander(pop_results)
write.csv(pop_results, file = here("D. documents", "appendix_table_c1.csv"))



########################## #
# Appendix D – Additional analyses of appropriateness ratings ####################
########################## #

# aggregate and reshape appropriateness ratings pooled by country
results_c <- df %>%
  mutate(condition = fct_recode(condition, Minority = "treatment",
                                Control = "control")) %>%
  group_by(condition, country) %>%
  summarize(
    vote_m = mean(app_vote, na.rm = T),
    vote_se = se(app_vote),
    consensus_m = mean(app_consensus, na.rm = T),
    consensus_se = se(app_consensus),
    leader_m = mean(app_leader, na.rm = T),
    leader_se = se(app_leader),
    chance_m = mean(app_chance, na.rm = T),
    chance_se = se(app_chance)
  ) %>%
  gather(variable, value, vote_m:chance_se) %>%
  separate(variable, c("variable", "stat"), sep = "_") %>%
  mutate(
    variable = fct_relevel(variable, "chance", "leader", "consensus"),
    country = fct_recode(
      country,
      DNK = "Denmark",
      HUN = "Hungary",
      IND = "India",
      RUS = "Russia"
    )
  ) %>%
  spread(stat, value)

# ~ Figure D1 ##########

# limit to vote and consensus appropriateness
ggplot(
  subset(results_c, variable == "consensus" | variable == "vote"),
  aes(
    x = fct_rev(condition),
    y = m,
    colour = fct_rev(variable)
  )
) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = m - se, ymax = m + se),
                position = position_dodge(width = 0.4),
                width = 0) +
  scale_colour_grey(start = 0,
                    end = 0.7,
                    guide_legend(title = "Decision rule:")) +
  xlab("") +
  ylab("Mean appropriateness ratings") +
  # ylim(c(-2.4, 2.4)) +
  facet_grid(country ~ .) +
  coord_flip() +
  theme_minimal() +
  theme(legend.position = "bottom")
ggsave(here("D. documents", "fig_d1.jpg"),
       width = 6,
       height = 4.4)


summary(aov(app_vote ~ country, df[df$condition == "control",]))

# calculate country level treatment effects

# changes in vote
by(df_pool, df_pool$country, function(x)
  country.t(x, app_vote))

# changes in consensus
by(df_pool, df_pool$country, function(x)
  country.t(x, app_consensus))

# ~ Figure D2 ##########

# aggregate and reshape appropriateness ratings split by country
results <- df %>%
  mutate(condition = fct_recode(condition, Minority = "treatment",
                                Control = "control")) %>%
  group_by(scenario, condition, country) %>%
  summarize(
    vote_m = mean(app_vote, na.rm = T),
    vote_se = se(app_vote),
    consensus_m = mean(app_consensus, na.rm = T),
    consensus_se = se(app_consensus),
    leader_m = mean(app_leader, na.rm = T),
    leader_se = se(app_leader),
    chance_m = mean(app_chance, na.rm = T),
    chance_se = se(app_chance)
  ) %>%
  gather(variable, value, 4:11) %>%
  separate(variable, c("variable", "stat"), sep = "_") %>%
  mutate(
    variable = fct_relevel(variable, "chance", "leader", "consensus"),
    country = fct_recode(
      country,
      DNK = "Denmark",
      HUN = "Hungary",
      IND = "India",
      RUS = "Russia"
    )
  ) %>%
  spread(stat, value)

# limit to vote and consensus appropriateness
ggplot(subset(results, variable == "consensus" | variable == "vote"),
  aes(x = fct_rev(condition), y = m, colour = fct_rev(variable))) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = m - se, ymax = m + se),
                position = position_dodge(width = 0.4),
                width = 0) +
  scale_colour_grey(start = 0,
                    end = 0.7,
                    guide_legend(title = "Decision rule:")) +
  xlab("") +
  ylab("Mean appropriateness ratings") +
  # ylim(c(-2.4, 2.4)) +
  facet_grid(country ~ scenario) +
  coord_flip() +
  theme_bw() +
  theme(legend.position = "bottom")
ggsave(here("D. documents", "fig_d2.jpg"),
       width = 6,
       height = 4.4)

# ~ Figure D3 ########
# limit to leader and chance appropriateness
ggplot(
  subset(results, variable == "leader" | variable == "chance"),
  aes(
    x = fct_rev(condition),
    y = m,
    colour = fct_rev(variable)
  )
) +
  geom_point(position = position_dodge(width = 0.4)) +
  geom_errorbar(aes(ymin = m - se, ymax = m + se),
                position = position_dodge(width = 0.4),
                width = 0) +
  scale_colour_grey(start = 0,
                    end = 0.7,
                    guide_legend(title = "Decision rule:")) +
  xlab("") +
  ylab("Mean appropriateness ratings") +
  # ylim(c(-2.4, 2.4)) +
  facet_grid(country ~ scenario) +
  coord_flip() +
  theme_bw() +
  theme(legend.position = "bottom")
ggsave(here("D. documents", "fig_d3.jpg") ,
       width = 6,
       height = 4.4)


# ~ Multilevel linear regression models for the appropriateness ###############

# baseline model with only random intercepts
mlm3.1 <- lmer(app_vote ~ 1 + condition +
                 (1 | ID) +
                 (1 | country_scenario),
               data = df)

# main model with random intercepts + random slope for treatment effect
mlm3.2 <- lmer(app_vote ~ 1 + condition +
                 (1 | ID) +
                 (1 + condition | country_scenario),
               data = df)
anova(mlm3.1, mlm3.2) # make sure model fit is improved
summary(mlm3.2)

# Regress Appropriateness of consensus

# baseline model with only random intercepts
mlm4.1 <- lmer(app_consensus ~ 1 + condition +
                 (1 | ID) +
                 (1 | country_scenario),
               data = df)

# main model with random intercepts + random slope for treatment effect
mlm4.2 <- lmer(app_consensus ~ 1 + condition +
                 (1 | ID) +
                 (1 + condition | country_scenario),
               data = df)
anova(mlm4.1, mlm4.2) # make sure model fit is improved
summary(mlm4.2)

# This code reproduces Table D1
# export results of models 3 and 4 in a nice table
stargazer(
  mlm3.2, mlm4.2,
  type = "text",
  covariate.labels = c("Intercept", "Vulnerable minority"),
  column.labels = c("(1)", "(2)"),
  model.numbers = F,
  intercept.bottom = F,
  ci = T,
  dep.var.labels = c("Vote", "Consensus"),
  title = "Table D1. Multilevel models for the appropriateness of voting and consensus",
  out = here("D. documents/appendix_table_d1.html")
)


######################################## #
### Appendix E – Main analysis without excluding comprehension failures ###################
######################################## #

nrow(df_all) # number of observations

# ~ Figure E1. #################

# calculate the proportion of each answer
df_prop_all <- df_all %>%
  filter(!is.na(decision)) %>%  # exclude NAs from decision variable
  mutate(condition = fct_recode(condition, Control = "control",
                                Minority = "treatment")) %>%
  group_by(condition, country) %>%
  summarize(
    vote = sum(decision == "vote"),
    consensus = sum(decision == "consensus"),
    leader = sum(decision == "leader"),
    chance = sum(decision == "chance"),
    trials = n()
  ) %>%
  gather(decision, frequency, vote:chance) %>%
  mutate(
    decision = fct_relevel(decision, "consensus", "leader", "chance"),
    country = fct_recode(
      country,
      DNK = "Denmark",
      HUN = "Hungary",
      IND = "India",
      RUS = "Russia"
    ),
    p = frequency / trials
  )


# This code reproduces Figure E1 in the appendix
ggplot(df_prop_all, aes(
  x = fct_rev(condition),
  y = p,
  fill = decision
)) +
  geom_col() +
  facet_grid(country ~ .) +
  # manual color codes help with unconventional order
  scale_fill_manual(values = c("#606060", "#bcbaba", "#969696", "#353535")) +
  scale_y_continuous(breaks = seq(0, 1, 0.25), labels = scales::percent) +
  coord_flip() +
  xlab("") +
  ylab("") +
  theme_minimal() +
  guides(fill = guide_legend(reverse = T)) +
  theme(
    legend.title = element_blank(),
    legend.position = "bottom",
    legend.spacing.x = unit(0.6, 'cm'),
    panel.grid.major.y = element_blank(),
    axis.text.x = element_text(size = 8)
  )
ggsave(here("D. documents", "fig_e1.jpg"), height = 3.8, width = 6.5)

# recalculate main aggregate changes in forced DV
df_pool_all <- df_all %>%
  group_by(ID) %>%
  mutate(
    vote = mean(decision == "vote", na.rm = T),
    consensus = mean(decision == "consensus", na.rm = T),
    leader = mean(decision == "leader", na.rm = T),
    chance = mean(decision == "chance", na.rm = T),
    app_vote = mean(app_vote, na.rm = T),
    app_consensus = mean(app_consensus, na.rm = T),
    app_leader = mean(app_leader, na.rm = T),
    app_chance = mean(app_chance, na.rm = T)
  ) %>%
  filter(!duplicated(ID))

# Recalculate the main results
# proportion of choices in control and treatment groups.
t.test(vote ~ condition, df_pool_all)
t.test(consensus ~ condition, df_pool_all)